home *** CD-ROM | disk | FTP | other *** search
/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc B) / Acorn User China CD-ROM (UK) (Disc B).bin / STUTTGART / LANG / GOFER / scripts_Power < prev    next >
Encoding:
Text File  |  1994-01-03  |  4.2 KB  |  130 lines
  1. -- Power series
  2.  
  3. data Power a = Series [a]
  4.  
  5. instance Add a => Add (Power a) where
  6.   (Series (x:xs)) + (Series (y:ys)) = Series ((x+y):zs) where
  7.        Series zs = (Series xs) + (Series ys)
  8.   (Series (x:xs)) - (Series (y:ys)) = Series ((x-y):zs) where
  9.        Series zs = (Series xs) - (Series ys)
  10.   negate (Series (x:xs)) = Series ((-x):ys) where 
  11.        Series ys = - (Series xs)
  12.   zero = Series zeros 
  13.  
  14. instance LeftMul a b => LeftMul a (Power b) where
  15.    a * (Series bs) = Series [ a*b | b <- bs ]
  16.  
  17. instance Div a b => Div (Power a) b where
  18.    (Series as) / b = Series [ a/b | a <- as ]
  19.  
  20. instance (LeftMul a b, Add b) => LeftMul (Power a) (Power b) where
  21.    (Series as) * (Series bs) = Series cs where
  22.       cs = [ sum [ a*b | (a,b) <- zip (rtake n as) bs ] | n <- [1..] ]
  23.  
  24. instance (Mult a, Add a) => Mult (Power a) where
  25.     unit = const unit
  26.  
  27. const :: Add a => a -> Power a
  28. const a = Series (a:zeros)
  29.  
  30. zeros :: Add a => [a]
  31. zeros = zero:zeros
  32.  
  33. vX :: (Mult a, Add a) => Power a
  34. vX = Series (zero:unit:zeros)
  35.  
  36. instance ( Div (Power a) a, LeftMul (Power a) (Power a),
  37.            Mult (Power a), Add (Power a), 
  38.            Add a) => Div (Power a) (Power a) where
  39.    xs / ys = xs * (invert ys)
  40.  
  41. invert :: (Mult (Power a), Add (Power a),
  42.            Div (Power a) a) => (Power a) -> (Power a)
  43. invert x = (Series (unit:ys))/x0
  44.     where  Series (xs@(x0:_)) = x
  45.            ys = [ f(i) | i <- [1..] ]
  46.            f(n) = -(xs!!n)/x0 - sum [f(n-i)*(xs!!i)/x0 | i <- [1..n-1]] 
  47.  
  48. {- The order of a nonzero power series is the least i for which the
  49.    i-th coefficient is nonzero. An infinite list of power series can be 
  50.    summed if the orders form a sequence that tends to infinity. From
  51.    the computational point of view we have to be given some function
  52.    g :: Int -> Int such that g(n) is a lower bound for the order of the
  53.    n-th power series. -}
  54.  
  55. limSum :: Add a => (Int -> Int) -> [Power a] -> Power a
  56. limSum g ps = Series [ f(n) | n <- [0..] ] where
  57.      f(n) = sum [ h n i | i <- [0..g(n)] ]
  58.      h n j = let Series zs = ps!!j in zs!!n
  59. -- i > g(n) && m < n ==> h m i == zero
  60.  
  61.  
  62. ptake :: Int -> (Power a) -> [a]
  63. ptake n (Series xs) = take n xs
  64.  
  65. -- take and reverse
  66. rtake :: Int -> [a] -> [a]
  67. rtake n xs = f n xs [] where
  68.    f 0 _ as          = as
  69.    f (n+1) (x:xs) as = f n xs (x:as)
  70.  
  71. zip :: [a] -> [b] -> [(a,b)]
  72. zip    []    _      = []
  73. zip     _    []     = []
  74. zip (a:as) (b:bs)   = (a,b):zip as bs
  75.  
  76. -- repeated differences
  77. diffs :: (Div (Power a) Int,
  78.           Add (Power a),
  79.           Mult ((Power a,Int) -> (Power a,Int))) => [a] -> [a]
  80. diffs xs = [ x | ( Series (x:_),_) <-
  81.                  [(next^n) (Series xs,0) | n <- [0..]]] 
  82.       where next (as,m) = let m' = m+1 in (((divX as) - as) / m',m')
  83.  
  84. divX :: (Power a) -> (Power a)
  85. divX   (Series as) = Series as' where _:as' = as
  86.  
  87. rtake0 :: Add a => Int -> [a] -> [a]
  88. rtake0 n xs = (rtake n xs) ++ zeros
  89.  
  90. undiff :: (LeftMul Int (Power a),
  91.            Add (Power a)) => Int -> [a] -> [a]
  92. undiff n xs = ys where 
  93.       Series ys = z - n*(timesX z)
  94.       z         = Series xs
  95.  
  96. timesX :: Add a => (Power a) -> (Power a)
  97. timesX (Series as) = Series (zero:as)
  98.  
  99. -- rebuild coefficients
  100. undiffs :: (LeftMul Int (Power a),
  101.            Add (Power a)) => Int -> Int -> [a] -> [a]
  102. undiffs d 0 xs     = rtake0 (d+1) xs
  103. undiffs d (n+1) xs = (take (n+2) (undiff (d-n-1) (undiffs d n xs)))
  104.                      ++ rtake0 (d-n-1) xs
  105.  
  106. --  interpolation for integral polynomials
  107. coeffs :: Int -> (Int -> Int) -> [Int]
  108. coeffs d f = rtake (d+1) cs where
  109.          cs = undiffs d d (diffs [ f(n) | n <- [0..] ])
  110.  
  111. print :: (Text a, Add a, Mult a, Ord a) => Char -> a -> [a] -> String
  112. print c x [] = show x
  113. print c x xs | x == zero = pr c 1 xs
  114.              | otherwise = (show x)++(spr c 1 xs)
  115.     where
  116.       pr c _ []                 = "0"
  117.       pr c n (x:xs) | x == zero = pr c (n+1) xs
  118.                     | otherwise = (t x c n)++(spr c (n+1) xs)
  119.       spr c _ []                = ""
  120.       spr c n (x:xs) | x == zero = spr c (n+1) xs
  121.                      | otherwise = (if x < zero then "" else "+")++
  122.                                    (t x c n)++(spr c (n+1) xs)
  123.       t x c n = (if x == unit then "" else
  124.                 if x == -unit then "-" else show x)++[c]++
  125.                 (if n == 1 then "" else "^"++show(n))
  126.  
  127. -- pretty print 
  128. display d f = print 'x' c0 cs where
  129.        c0:cs = coeffs d f
  130.